Reliability of temporal coding on pulse-coupled networks of oscillators 
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We study the reliability of spike output in a general class of pulse-coupled oscillators receiving a 
fluctuating input. Showing that this problem is equivalent to noise-induced synchronization between 
identical networks of oscillators, we employ the phase reduction method to analytically derive the 
average Lyapunov exponent of the synchronized state. We show that a transition occurs between 
reliable and unreliable responses at a critical coupling strength, which is determined through the 
competition between the external input and recurrent input. To our surprise, the critical value does 
not depend on intrinsic properties of oscillators. 
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Noise-induced synchronization appears in a variety of 
phenomena including lasers chemical reactions 0, 
gene networks [Sj] and neuronal systems 0, Q . In these 
systems, periodic or chaotic oscillators driven by a com- 
mon fluctuating input synchronize with each other due to 
the nonlinearity of oscillators and the stochastic nature of 
the input Q . The phase reduction Q and the Lyapunov 
analysis proved that two or more identical oscillators re- 
ceiving a common fluctuating input are always in-phase 
synchronized regardless of their intrinsic properties and 
initial phases [I,!!] . We can interpret such oscillators as a 
single oscillator receiving the same input repeatedly, but 
with different initial phases, i.e. many trials of an input 
application. Therefore, the in-phase synchronization of 
input-driven oscillators implies, in a single oscillator, the 
reproducibility of the responses to a repeated input, or 
response reliability, which is particularly important for 
processing external signals. Reliable responses to a fluc- 
tuating input are actually measured from single cortical 
neurons \4^. However, neurons and other oscillators in 
the real world work collectively in their networks rather 
than individually. To study whether a network of os- 
cillators still has response reliability, we develop a the- 
ory of noise-induced synchronization between networks 
of oscillators rather than between single oscillators. We 
find a transition from reliable to unreliable responses at 
a critical coupling strength. Deriving average Lyapunov 
exponent analytically, we reveal that the critical value 
is determined through the competition between variance 
of the external input and of internal recurrent inputs re- 
gardless of details of oscillators. Around the transition 
point where magnitude of the average Lyapunov expo- 
nent is small, information of initial states can stay in the 
network for a long time. We discuss a possible role of the 
long time scale in role-sharing between rate and temporal 
coding on neuronal computation in the brain. 



A network of pulse-coupled N limit-cycle oscillators 



receiving fluctuating inputs are described as: 

N 



dt 



F{X,) + Ut) + E E a^,5{t ~ tf^^^)g, (1) 

J — 1 spike 



where i — 1, - ■ ■ ,N and X = F (X) has a stable limit- 
cycle solution Xq (t) . A unit vector g indicates the 
direction of interactions in the multidimensional space 
spanned by X. We assume, with neuronal oscillators 
in mind, coupling matrix gij is a sparse random matrix 
with connection probability p, and each nonzero com- 
ponent of gij is either g or —gj if j refers to an exci- 
tatory or an inhibitory cell, respectively. The network 
consists of Ne excitatory neurons and Nj — N — Ne 
inhibitory neurons. Fluctuating external inputs repre- 
sent independent white Gaussian processes with strength 
(li(t)lj(s)) — 2a'^dij6 {t ~ s), and tj^'^'^^ represents 
spike times of the jth neuron. We use g as a con- 
trol parameter of the network and require, for simplic- 
ity, that gi is in proportional to g and satisfies the bal- 
ance condition [10, gNs + {—gi)Nj = 0, whereas 
results of the paper are independent of the restriction. 
When g — 0, response of oscillators are always reli- 
able, i.e. the spike sequence of each oscillator con- 
verges into the same sequence in different trials. Figure 
[T^ and [TJj demonstrates the reliable responses obtained 
from numerical calculations of quadratic integrate-and- 
fire (QIF) neurons, F (X) = + /, with variable re- 
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Whereas two trials start from different initial states, 
raster plots of them converge into same sequences, i.e. 
same spike times. We then introduce finite couplings 
(7 > to the network of oscillators and calculate firing 
responses in a similar way to Fig. [T|d. When coupling 
strength is small, the population is still reliable, spike 
sequences of different trials converge into the same one 
(Fig. [1];). However, the reliability is lost from the pop- 
ulation when coupling strength is sufficiently large (Fig. 
[TJi). Spike sequences of two trials never converge into 
the same one while we apply the same input to trials. In 
terms of synchronization, the result means that fluctuat- 
ing inputs induce phase synchronization to two identical 
networks of oscillators only when coupling strength of the 
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FIG. 1: Responses reliability of a pulse-coupled network of 
quadratic integrate-and-fire models. I = 1.0, a — 1.0, Ne ~ 
80 and — 20 (a) Time evolutions of X of an oscillator for 
1st (black) and 2nd (gray) trials when g = 0. (b-d) Raster 
plots of spike times of randomly chosen 10 oscillators of the 
network for 1st (black) and 2nd trials (gray), g = (b), 0.1 
(c) and 0.3 (d). 



network is sufficiently weak. 

Regarding the fluctuating signals and recurrent con- 
nections as perturbations to the deterministic oscillators, 
we apply the standard phase reduction method to Eq. ([1]) 
and obtain stochastic equations of phases (f>i as. 



N 

j — 1 spike 



where w is an intrinsic frequency of the unperturbed os- 
cillators. Phase sensitivity function or phase response 
function is defined uniquely from Xo{(f>) as Z{(f)) — 
gradx0|x=Xo(0) ,''']■ To simplify notations, we assumed 
without loss of generality that g is in parallel with ^ and 
replaced vector variables Z (0) g and Z (0) ^ to scalar 
variables Z (0) and Z {(j)) ^ in Eq. ©. For the QIF 
model, for instance, w = /~^/^, Z {(j)) = uj {1 + cos (0)). 
Note that discontinuous variable resetting X {^t^v^'"^^ = 
—00 is now reduced to continuous 
dynamics over a spike threshold </> = tt because deriva- 
tion of F is continuous between X = 00 and X = — 00, 
F'(oo) — F'{—oo). To avoid unrealistic cases where an 
oscillator sends numerous spikes within a short interval of 
time when its fluctuating phase crosses the firing thresh- 
old, Z {(j)) should be vanish around the threshold. Real- 
istic neuron models including the QIF model satisfy the 
condition. Reliability of the firing responses is equivalent 
to that of phase dynamics because phase deviation is pro- 
portional to deviation of firing time. Phases of two trials, 
and evolve satisfying Eq. ^ from different initial 
phases but receiving same inputs, ^i, in the same net- 



work, gij. Since phase synchronized state, (pi = 4>i^ is an 
obvious solution of these two equations, linear stability 
around the solution determines the response reliability. 
To evaluate the stability we linearize Eq. ^ in terms of 
small phase differences, ipi = 4'i ~ and calculate av- 
erage Lyapunov exponent [l^ over all oscillators in the 
network. 

Coupling terms are linearized as follows. Consider in- 
crement of ipi, when jth cell fires. We can take firing time 
of ist = 0, firing time of 0^ is therefore t — dT — ipj/^- 
We can assume dT > without loss of generality. At 
t = 0, 4>i receives a spike from 0^, 



0,(0+) = M0) + 9^JZiM0)) 

0,(0+) = 0,;(O) 



(3) 



Phases evolve as follows from t = 0+ to dT, because dT 
is a short interval, 

0,(dr) - 0,(0+) + {uj + a2Z'Z(0,(O+))) dT 

+Z{UO+))dW , (4) 

0,(dr) = 0,(O+) + (u; + a2Z'Z(0,(O+)))dT 
+Z{UO+))dW 

where dW — ^idT. In order to evaluate phase responses 
Z(0) at precise timings just before spike inputs, we trans- 
lated Eq. (12) to equivalent Ito integrals in Eq. (g]) p^ . 
Third terms of Eq. ^ result from the translation. Fi- 
nally, at t — dT, (pi receives a spike from 0j, 



0,(dT+) = 0,(rfT) 

0,(dT+) = 0,(dT)+gyZ(0,(dT)) 



(5) 



^_sp^ke S| I ^ -y^g linearize from Eq. ([3]) to (O in terms of ip with 

an attention that dW is the order of dT^^'^ and then the 
1/2 . 

(2) order of tpj ■ Taking all connections into account and 



neglecting terms higher than the order of g^, we obtain 
linearized equation of ip as 



dipi 
dt 



(Z'(0,)Z(0,))' + Z'(0,)e,U, 



(6) 



N 



j—1 spike 

By introducing new variables iji = {logipf) /2, Eq. ^ is 
further rewritten as 



^ = a\Z'{^,)Z{<t>i))' + Z'{<Pi)^, 



o^Z'i 



(7) 



N 



J — 1 spike 



l + gyZ'(0,) 1 



ipi 



Since the Lyapunov exponent Ai is defined as 
\miT^oo{Vi{T) — yi{Oj)/T, the long time average of the 
Eq. ([7]) coincides with A^. We assume that the network 
is in asynchronous steady firing state due to fluctuating 
inputs and replace spike times of cells by independent 
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Poisson processes with firing rate r [15|. Then averag- 
ing of Eq. ([7|) over the Poisson processes and over all 
oscillators in the network gives 

where (Z'^) = {2tt)-^ J^"" Z' {(t)f dc/). Here we used the 
assumption of weak inputs and weak interactions and re- 
duced distributions of phases to uniform distributions in 
[0, 2tt]. Unfortunately, Eq. ^ is not a closed form of y^. 
However, when variance of ip"^ is small, or ~ (V'^)^i 
the last term of Eq. ^ is approximated as {tpf /^p-j) = 1 
and we finally obtain the following main formula of the 
average Lyapunov exponent: 

A = i-a' + Nr{g')) . (9) 

Note that when variance of -0^ is not small, Eq. ([9]) gives 
the lower bound of A because {ipf/i^'j) > 1 in generally. 

To confirm the above analysis, we calculate averaged 
dynamics of yi = (log?/'^) /2 numerically for networks of 
QIF oscillators. Due to fluctuating inputs and recurrent 
interactions, yi themselves do not evolve monotonically. 
However, population averages of yi decrease or increase 
almost linearly depending on coupling strengths g as pre- 
dicted by Eq. dH). 

Our expression of the Lyapunov exponent, Eq. 
tells us two important facts of the reliability transition. 
First, the transition stems from a competition between 
two variances, variance of input signals and variance of 
recurrent inputs Nr{g'^). Whereas the first contribution 
to the exponent is negative, the second is always posi- 
tive. Therefore, the network lost their reliability when 
the second exceeds the first. Second, the factor of {Z'"^) 
which reflects intrinsic properties of oscillators is mul- 
tiplied equivalently to these two factors and Nr{g'^) 
in Eq. ([9]). Therefore, the critical coupling strength gc, 
given as the solution of A = 0, is universal in the sense 
that gc is independent of details of oscillators. For the 
network structure we used in Fig. [1] the critical value is 
given as gc = (J \J {Npr)~^Ni/NE regardless of oscillators 
on the network. If we use another natural normalization 
of coupling strengths as g ^ g/^/Np 11], we can elimi- 
nate N from the critical strength, gc = a^Jr~'^Ni /Ne- 

In the vicinity of the critical coupling strength g = gc 
where |A| <C 1, information on the initial states of oscilla- 
tors may disappear quite slowly after the onset of input. 
This slow transient behavior might have the following 
implications for computations by cortical networks. The 
output of the computation is not simply determined by 
the current input, but is also modulated by the brain's 
internal state and/or input histories [l6l . Il7l |. In our 
model, the membrane time constant sets the short time 
scale that enables the network to respond quickly to an 
external input with firing rate of population dynamics 
[ill. ITsl. By contrast, the critical dynamics of temporal 
spike sequences may set a much longer time scale to en- 
sure the response diversity reflecting the initial state or 
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FIG. 2: Time evolutions of logarithm of phase differences yi 
of randomly chosen 5 oscillators (gray lines). Sudden tenta- 
tive jumps to large values on these lines are because of differ- 
ences of spike times between two trials (see Eq. ([Sjl and ((Sjl). 
Thick black curves are averages of yi over all oscillators in the 
network. Slopes of these curves agree well with the analytical 
results of Eq. ((9)l shown by slopes of thin black lines. The 
values of parameters are the same as in Fig. [T] except that 
a = 0.2, and g = 0.01 (a), 0.05 (b) and 0.1 (c). 



input histories. This implies that neuronal populations 
may simultaneously achieve two different time scales by 
parallel use of rate code and temporal code. Further 
studies are required for clarifying this possibility. 

So far, we have restricted our study to super-threshold 
neurons which continue to fire without external inputs. 
Numerical simulations of sub-threshold QIF model with 
/ < 0, however, suggest that similar transition also oc- 
curs in a network of sub-threshold neuron models (fig- 
ure [3]). It remains unknown whether this transition may 
appear in a broad class of sub-threshold neuron mod- 
els because the phase reduction method is not applicable 
to sub-threshold neuron models. A unified treatment of 
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FIG. 3: Responses reliability of a network of sub-threshold 
QIF models for two trials. The values of parameters are the 
same as in Fig. [T]except that I — —1.0, and g = 0.4 (a) and 
1.0 (b). 

super- and sub- threshold neuron models is awaited. We 
have assumed that couplings among oscillators are delta- 
functions. To remove a doubt that our results might be 
pathological phenomena come from singularity of delta- 
functions, we calculated reliability of coupled oscillators 
numerically using a alpha-function, a(t) = a^i exp(— ai), 
instead of 5{t). Again, we could see similar transition 
from reliable to unreliable responses (results not shown) . 
Linear integrate-and-fire model is the most useful de- 
scription of firing neurons. This model, however, behave 
unrealistically about response reliability even when g = 
because of its anomalous variable resetting Here, 
we use quadratic integrate-and-fire model to avoid the 



problem. Coupled oscillators may synchronize with each 
other if g is sufficiently large, whereas we have only con- 
centrated on the asynchronous steady state. Synchro- 
nization may affect average Lyapunov exponent and may 
change the transition significantly because we must use 
correlated stochastic processes instead of independent 
Poisson processes when we average Eq. (jS]) to obtain 
A. 

In conclusion, coupled reliable elements are not nec- 
essarily reliable any more. Spike responses of coupled 
oscillators to fluctuating inputs show transition from reli- 
able responses to unreliable responses. In terms of noise- 
induced synchronization, common noises fail to induce 
phase synchronization to networks of strongly coupled 
oscillators whereas same inputs always induce synchro- 
nization to single oscillators. Underlying mechanism of 
the transition is competition between a variance of ex- 
ternal signals and a variance of internal recurrent inputs. 
Critical coupling strength derived analytically is inde- 
pendent of details of oscillators because phase response 
functions appear equivalently in these competing factors. 
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